Worm-type Monte Carlo simulation of the Ashkin- Teller model on the triangular 

lattice 
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We investigate the symmetric Ashkin- Teller (AT) model on the triangular lattice in the antifer- 
romagnetic two-spin coupling region (J < 0). In the J -co limit, we map the AT model onto a 
fully-packed loop-dimer model on the honeycomb lattice. On the basis of this exact transformation 
and the low-temperature expansion, we formulate a variant of worm-type algorithms for the AT 
model, which significantly suppress the critical slowing-down. We analyze the Monte Carlo data by 
finite-size scaling, and locate a line of critical points of the Ising universality class in the region J < 
and K > 0, with K the four-spin interaction. Further, we find that, in the J — > — oo limit, the critical 
line terminates at the decoupled point K = 0. From the numerical results and the exact mapping, 
we conjecture that this 'tricritical' point (J — > —oo,K = 0) is Berezinsky-Kosterlitz-Thouless-like 
and the logarithmic correction is absent. The dynamic critical exponent of the worm algorithm is 
estimated as z = 0.28(1) near (J — > -co, K — 0). 

PACS numbers: 



I. INTRODUCTION 

The Ashkin- Teller (AT) model is a generalization of 
the Ising model to a four-component system of which 
each lattice site is occupied by one of the four states [J- 
7]. In 1972, Fan [l| associated each lattice site with two 
Ising variables (er, r) and represented the four states by 
the combined states (1, 1), (1, —1), (—1, 1) and (—1, —1). 
On this basis, the reduced Hamiltonian (fc^T = 1) of the 
AT model reads 
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where the sum (ij) runs over all the nearest-neighbor 
pairs of spins, J CT ( J T ) represents the two-spin interaction 
for a (t), and K is the four-spin interaction. Examples of 
physical realizations of the AT model include: 1) , systems 
with layers of atoms and molecules adsorbed on clean 
surfaces-e.g., selenium adsorbed on the Ni(100) surface 
[U and oxygen-on- graphite system and 2), systems 
with layers of oxygen atoms in the CuO plane, like high- 
T c cuprate YBCO 

The AT model exhibits very rich critical behavior and 
plays an important role in the field of critical phenomena. 
Figure Q] displays the phase diagram of the AT model on 
the square lattice for J = J a = J T > (we shall only 
consider this symmetric case in this work). The model 
reduces to two decoupled Ising systems for K — 0, and 
is equivalent to the 4-state Potts model along the diago- 
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FIG. 1: (Color online) Phase diagram of the AT model on the 
square lattice. The 'P-I-O' curve (thick cyan line) is self-dual 
and has continuously varying critical exponents, separating 
the paramagnetic and the ferromagnetic state in Ising vari- 
ables a, r, and ar. The 'P-A', 'P-B' and 'O-C lines are 
commonly believed to be Ising-like, which are represented by 
thin magenta lines. 



nal line J ~ K. The whole 'P-I-O' line is critical, with 
continuously varying critical exponents, and with the de- 
coupled Ising point I and the 4-state Potts point P as 
two special points. The two branches 'P-A' and 'P-B' 
are also critical, and are numerically shown to be in the 
Ising universality class. On other two-dimensional pla- 
nar lattices like the honeycomb, triangular, and kagome 
lattices, the phase diagram of the AT model with J > 
is similar as Fig. [TJ except the fact that the antifcrro- 
magnetic transition line for K < may be absent on 
non-bipartite lattices like the triangular and the kagome 
lattice. 
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In this work, we shall consider the AT model on the 
triangular lattice. From the duality relation and the star- 
triangle transformation, it was already found [ll| in 1979 
that the critical P-I-0 line is described by 

e- iK = \(e iJ -l), (2) 

with K < i log 2. Further, it can be shown that the 
model on the infinite-coupling point 0(J = — K — > oo) 
can be mapped to the critical O(n) loop model with n = 2 
on the honeycomb lattice, and the well-known Baxter- Wu 
model on the triangular lattice at criticality [l2| . In the 
limit J = K — > — oo, the model is equivalent to the 4- 
state Potts antiferromagnet at zero temperature, which is 
also critical. In the limit J = 0, K — > — oo, the AT model 
reduces to the zero-temperature Ising antiferromagnet in 
variable or; the same applies to the limit K — 0, J — > 
— oo for the two decoupled Ising variables a and r. Phase 
transition of the triangular-lattice Ising antiferromagnet 
is absent at finite temperature, and at zero temperature 
the system has non-zero entropy per site [H, The 
pair correlation on any of the three sublattices of the 
triangular lattice decays algebraically as a function of 
distance, and the associated magnetic scaling dimension 
isX h = l/4 [H. 

On the square lattice, the phase diagram for J < is 
the symmetric image of Fig. [T] with respect to the K axis 
(J — > — J), arising from the bipartite property. However, 
to our knowledge, the phase diagram of the AT model 
is still unknown on the triangular and the kagome lat- 
tice with J < 0. Clearly, the Ising critical line P — A 
should continue into the region J < 0,K > 0, albeit 
it remains to be explored how this extension looks like. 
Due to the absence of exact result, we will apply Monte 
Carlo method and the finite-size scaling theory. Monte 
Carlo simulation of the triangular AT model is challeng- 
ing for large negative coupling J < 0, arising from the so- 
called geometric frustration. Antiferromagnetic coupling 
J < means that the neighboring Ising spins prefer to 
be anti-parallel. However, such a preference cannot be 
satisfied for all of the three neighboring pairs on any ele- 
mentary triangular face. One can at most have two anti- 
ferromagnetic pairs. For such a frustrated system, most 
Monte Carlo simulation suffers significantly from criti- 
cal slowing-down. In fact, as J — » — oo, the Metropolis 
and the Swendsen- Wa ng-type cluster algorithm are found 
to be non-ergodic 0, [IoHl8l |. Recently, worm-type al- 
gorithms with the so-called rejection- free was developed 
for the antiferromagnetic I sing model on the triangular 
lattice and other systems [l9j, |20|. This algorithm has 
been proved to be ergodic at zero temperature and only 
suffers from minor critical slowing-down. The rejection- 
free worm algorithm can be extended to the AT model, 
albeit the efficiency is limited for nonzero K in the zero- 
temperature limit J — > — oo. 

The outline of this paper is as follows. Section [XT] de- 
scribes the partition sum of the AT model as well as 
an exact mapping to a fully-packed loop-dimer (FPLD) 
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FIG. 2: (Color online) A spin configuration of the AT model 
on the triangular lattice and the corresponding LT-expansion 
graph on the honeycomb lattice. Thick blue line represents 
blue bond; thin red line denotes red bond; the same below. 

model in J — > — oo limit. A variant of worm-type algo- 
rithms is developed in Sec. IIIII The numerical results 
are presented in Sec. [IVJ In Sec. [V] we investigate the 
dynamic critical behavior of one of the worm algorithms. 
A discussion is given in Sec. IVI[ including the phase di- 
agram on the kagome lattice. 

II. MODEL AND EXACT MAPPING 

A. Low-temperature expansion of the AT model 

Instead of directly updating the spins, the worm-type 
algorithms [2l|, [22[ for the Ising model simulate the 
graphical representation which can be the high- and the 
low-temperature expansion graphs. The worm methods 
in Refs. [U [22j can be generalized to the graphical ex- 
pansion of the AT model. In the following, we shall use 
the low-temperature (LT) expansion, defined on the dual 
lattice of the triangular lattice — the honeycomb lattice. 
Given a spin configuration {a, r}, for each pair of nearest- 
neighboring vertices one places on its dual edge: 

• nothing if <7j = <jj , = Tj , 

• a red occupied bond if Oi — Uj , Tj 7^ Tj , 

• a blue occupied bond if <Ji ^ <7j , Tj = Tj , 

• a red and a blue bond if <Ji ^ Oj , 7^ Tj . 

In other words, depending on the associated pair of spins 
on the triangular lattice, an edge on the honeycomb lat- 
tice can be in one of the four states: vacant, red, blue, 
and red+blue. An example is shown in Fig. [21 Since the 
coordination number is 3 for the honeycomb lattice, the 
red and blue bonds form a series of disjointed loops in 
red and blue color, respectively. Note that the red and 
the blue loops are allowed to share common edges. In 
this way, a spin configuration on the triangular lattice 
is mapped onto a loop configuration on the honeycomb 
lattice, while a loop configuration corresponds to 4 spin 
configurations [3l|, which are related to each other by 
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FIG. 3: (Color online) Vertex states in the LT-expansion 
graph of the AT model. 

globally nipping the a or/and r Ising spins. Let \E r \, 
\E{,\, and |£y+&| be the number of red, blue, and red+blue 
bonds, the partition sum of the AT model can be written 
as (up to an unimportant factor) 

Z A T = J2 X l ErlX l Eblx l E + l +bl , (3) 

{C} 

where the summation {£} is over all loop configurations. 
From the mapping, one can obtain the relative statistical 
weights as 

X r = X b = e - 2J - 2K and X r+b = e~ 4J . (4) 

One can further describe the AT model in the language 
of the vertex states, which will serve as the basis for the 
formulation of the worm-type algorithms in this work. 
In the loop configurations, all the vertices must have an 
even number of incident red (blue) bonds. Accordingly, 
only the 5 types of vertex states in Fig. [3] exist, where 
the states are unchanged under spatial rotations. Simple 
calculations yield the statistical weights as 

Wi = i , w 2 = w 3 = e - 2J - 2K , 

Wi = e~ 4J , and W 5 = e ~ AJ ~ 2K . (5) 

Let \Vi\ be the number of vertices at state i with i = 
1,2,3,4,5, the partition sum of the AT model can be 
written as (up to a constant) 

ZAT=j:f[wr i , (6) 

{V} i=l 

where the summation {V} is over configurations with ver- 
tex states in Fig. [3] 

B. Exact mapping in the J — ^ — oo limit 

Given a finite four-spin coupling K, when the antifcr- 
romagnetic coupling J becomes stronger and stronger, 
more and more vertices will be at state-4 and -5 in Fig. |3l 
because W± ~ W5 oc exp(— 4J) increases faster than 
W\,W%,Wz, as seen from Eq. ([5]). In the J — > — 00 limit, 
only state-4 and -5 survive. We can then redefine the 
edge states in state-4 and -5 as following. The empty 
edge is replaced by a 'dimer', while the 'blue+red' edge 
is regarded as empty; namely, the edge is now at state: 
empty, dimer, red, or blue. As a result, state-4 and -5 
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FIG. 4: (Color online) (a), State-4 and -5 after the redefinition 
of the edge states, (b), Vertex states in the FPLD model. The 
dashed black line represents dimer. 

become those in Fig. |U[a) . One observes that the occu- 
pied bonds at state-5 form a series of disjointed loops; 
these loops are now constructed by bonds alternatively 
in color red and blue. Further, one notes that the color- 
degree freedom can be simply integrated out, and each 
loop gains a statistical- weight factor 2. Without the color 
information, the edge is at state: empty, dimer, or bond, 
and the vertex states reduce to those in Fig. H^b), where 
new labels '6' and '7' are used. The statistical weights 
are 

W 6 = 1 , W 7 = e~ 2K . (7) 

On this basis, the partition sum of the AT model in 
the J — > — 00 limit can be written as 

Z FPLD =^r/wf 71 , (n = 2) (8) 

{V} 

where the summation is over configurations with all ver- 
tex states in Fig. |4j and £ is the number of loops. We 
shall refer to the model defined by Eq. (j8]) and Fig. [4] as 
the n-color FPLD model. 

Note that, for finite K, the loops in the FPLD model 
are 'dilute' due to the presence of state-6. However, in 
the K — > —00 limit, only state-7 survives, and one ob- 
tains the mapping between the triangular 4-state antifcr- 
romagnet at zero temperature and the honeycomb n = 2 
fully-packed loop model. For K — > 00, the model reduces 
to the fully-packed dimer model, which is equivalent to 
the triangular Ising antiferromagnet at zero temperature. 

We conclude this subsection by mentioning that the 
FPLD model is very similar to the honeycomb O(n) loop 
model 23]. The difference is that in the former the ver- 
tices off the loops are paired up by dimers, while not in 
the latter. Namely, the configuration space for the FPLD 
model is a subspace in the 0(n) loop model. Albeit it 
remains to be explored whether or not the two models 
are in the same universality class, it is not surprising if 
this turns out to be the case. 



III. WORM ALGORITHMS 

The worm algorithm for the high-temperature expan- 
sion graphs of the Ising model was first formulated by 
Prokof'ev and Svistunov [2l| . and the dynamic criti- 
cal behavior was studied in Ref. [22j. Recently, Wolff 
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FIG. 5: (Color online) Additional vertex states in the worm 
method for the AT model. The red (blue) filled circle denotes 
a defect in red (blue) vertex configuration. 



provided a worm-type simulation strategy for O(N) 
sigma/loop models 24]. The underlying physical pic- 
ture of the worm method is beautifully simple: enlarge 
the state space of the to-be-simulated model, define an 
extended model, and simulate the system by a local al- 
gorithm. 



A. Worm algorithm for finite J 

Let us now generalize the worm method in Refs. [2lL 
[22l | to the AT model in the language of the vertex states, 
defined by Eq. © and Fig. [5] 

Enlarge the state space. We first introduce new ver- 
tex states by deleting from (or adding to) the states in 
Fig. [3] a red or blue bond. This leads to the 8 additional 
vertex states in Fig. [5] The state space is then enlarged 
such that a configuration has a pair or none of vertices 
at states in Fig. [SJ Such a pair of vertices are named 'de- 
fects' and denoted as (u,v). Accordingly, the state space 
can be divided into two subspaces: one without defect 
(u = v) and the other with two defects («/»); we shall 
refer to them the M (measuring) and W (worm) sector, 
respectively. A careful check yields that the pair of de- 
fects in the W sector must be connected via a string of red 
or blue occupied bonds. Namely, u and v are either both 
at states {1', 4', 5', 7'} or {2', 3', 6', 8'} in Fig. GO For the 
later convenience, we let u, v be ordered as u <— v, and 
thus the interchange (u ■<-> v) would lead to a different 
configuration for u 7^ v . 

Define the extended model. With the inclusion of the 
defects and the vertex states in the W sector, a configura- 
tion can now be completely specified by its vertex states 
{V}, and the ordered pair of defects (u,v). The partition 
sum of the extended model can be separated into two 
parts. The part in the M sector is defined as 



1 5 

z M = z AT = — 6 U=V JJ W i 



(9) 



{V,u,v} 



where V is the volume of the system and 6 is the Kro- 
necker delta function. The summation {V} is over vertex- 



state configurations with states in Fig. [3] and coordina- 
tions (u,v). Factor 1/V accounts for the summation of 
u = v over the whole lattice. Similarly, the part of the 
partition sum in the W sector can be defined by 



■■w 



1 

V 



{V,u,v} 



(10) 



where the summation {V} is over configurations with two 
vertex states in Fig.[5]and all other in Fig.[3j and Wj> are 
the statistical weights for states in Fig. [5] The extended 
model can then be defined as 



fill 



with (; w > a constant factor controlling the relative 
weight in the M and the W sector. 

For a complete definition of the extended model, the 
statistical weights Wy for states in Fig. [5] should have 
a definite value. It is natural that they arc defined in 
accordance with the edge states, which lead to 



Wv =W 2 > = e 



-J-K 



W y = W v = e 



-3J-K 



W w = W& = e - 3J - 3K , W r = W 8 > = e~ 5J - K 



(12) 



Formulate the worm algorithm. One can now use any 
valid algorithm to simulate the model defined by Eq. (jTTjl . 
Since a configuration is specified by the ordered triplet of 
parameters (V, u, v), an update can be acted on the ver- 
tex states V and/or the locations of defects (u,v). The 
worm strategy is to randomly move u and/or v around 
the lattice and update V by changing the edge states dur- 
ing the biased random walk. Suppose that ti/« are in 
red (in the W sector). As u moves to a neighboring vertex 
u n , the edge (uu n ) state will be symmetrically updated: 
a red bond is placed (deleted) if it is absent (present). In 
this way, state at u will be back in Fig.[3]after u moves u n . 
Accordingly, the number of defects remains unchanged if 
v =/= u n or becomes zero if v = u n . This accounts for 
a step of random walk in the W sector or from the W 
to the M sector. For the case u = v, by the symmetric 
update of edge state, one will generate a pair of defects 
which can be either in red or blue. Therefore, one never 
introduces more than two defects. The parameter f = 1 
is set in this work, and a version of the worm algorithm 
reads (Algorithm 1) 

1. If u = v, randomly choose a new vertex vl and set 
u = v = v! . Equally choose color red or blue for 
the to-be-proposed defects; say red. 

2. Interchange u O v with probability 1/2. 

3. Randomly choose one neighboring vertex u n of u. 
Propose to move u — > u n . 

4. Propose to symmetrically update the edge-uu„ 
state: red ■(-> vacant and blue -H- red+blue. 
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5. Accept the proposal with probability 



V a = min 



h(Wi a) Wi a J)/(W^W^) 



according to the Metropolis-Hasting scheme. The 
superscript (b) and (a) means "before" and "after 
update", respectively. The statistical weights are 
given Eqs. © and ([12]) . 



Monte Carlo simulation of the AT model consists of repe- 
tition of these steps. The detailed balance at each step is 
straightforward since the algorithm is just a Metropolis- 
type update. If one regards the connected pair of 'de- 
fects' as a worm, the above steps mimic the crawling of 
the worm on the lattice. This is responsible for the ter- 
minology 'worm'. 

Measurement. Measurement can take place either in 
the whole enlarged state space or in the M subspace. For 
the high-temperature graph of the Ising model, it can 
be shown that the partition sum of the extended model 
is related to the Ising model as Z wmm — xZi S i ng , where 
X is the magnetic susceptibility. Thermodynamic quanti- 
ties can be measured in the enlarged configuration space. 
Nevertheless, if one is only interested in the original sys- 
tem, it is sufficient to sample in the M sector. This would 
define a Markov subchain with a coarse unit of Monte 
Carlo step between two subsequent configurations in the 
M sector. The detailed balance is clear since it is satisfied 
in each basic step in Algorithm 1. 

Improved version. As mentioned earlier, state-4 and 
-5 (Fig. \S§ would dominate in the M sector as J — > — oo; 
analogously, only state- 7' and -8' (Fig. [5]) survive in the 
limit, as seen from Eq. (|12[) . This implies that, as soon as 
both u and v are at state- 7' and —8', they will be frozen 
there forever, and thus Algorithm 1 becomes non-ergodic. 

The same difficulty occurs for the worm simulation 
of the triangular Ising antiferromagnet at zero temper- 
ature. A rejection-free technique was introduced 0, [2(| 
to overcome such a problem, based on the observation 
that the detailed balance in the coarse step does not re- 
quire the detailed balance in each basic step in the W 
sector. Let u n (n = 1, 2, 3) denote the neighbouring ver- 
tices of u and p n be the probability that u moves to u n 
in Algorithm 1, the probability for u to be unmoved is 
Po = 1 — (P1+P2+P3). The absorbing problem of (u, v) at 
state- 7' and -8' is reflected by po — > 1 as J —> —00. In the 
W sector, one can explicitly set zero for the probability 
that u remains unmoved, and defines the new transition 
probabilities p' n as 



Po 



Pl = P2 = P3 
Pi P2 P3 ' 

l-(p[+p' 2 +p' 3 ) = 



(13) 



The details can be found in Refs. |19l |20|. 

The absorbing problem can also be solved in the 
present formulation of the worm algorithm. Actually, 
the absorbing problem is somewhat 'artificial' here, since 
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FIG. 6: (Color online) Vertex states in the W sector for the 
FPLD model. The black filled circle denotes a defect. 



it arises from the particular assignment of the statisti- 
cal weights to states in Fig. [5] by Eq. (TT2"j) . There is no 
reason, however, why one should use Eq. (fT2l) if only the 
original AT model (0) is of interest. The absorbing prob- 
lem simply dissolves if the statistical weights are given 
by 

Wv = W 2 > = e - 2J - 2K , Wy = W v = e~ 4J , 

W 5 > = W 6 > = e - 2J - 2K , W v = W 8 , = e- 4J - 2K (14) 

Other definitions are possible. 



B. Worm algorithm for 

Algorithm 1 using Eq. (fT4")) is found to be efficient in 
most of the region with > J > — 00 and for small 
K in the J — > —00 limit. In this limit, the efficiency 
significantly drops as K deviates from 0. 

Hereby we shall make use of the exact mapping of the 
AT model onto the FPLD model © and formulate an- 
other version of the worm algorithm. Following the same 
procedure in the above subsection, we first introduce 5 
additional states in Fig. [BJ The partition sum in the M 
sector is defined as 



Z 



M 



Z 



FPLD 



1 

V 



E ^ 

{V,u,v} 



,,nW', 



V7I 



(15) 



with n = 2. Again, the summation is over configurations 
with states in Fig. @] and over the location of u = v. The 
partition sum in the W sector is given by 



- E 

V ^ 

{V,u,v} 



13 



(16) 



The extended model is defined by Eq. (ITT1) . 

The formulation of the worm algorithm follows the 
standard strategy in the above subsection, except that 
the edge-state update should take a different scheme. Let 
e = 0,1,2 denote the edge-e state 'empty', 'bond', and 
'dimer', respectively, and define the module-3 summa- 
tion rule as mod 3(e + Ae) with Ae = 1, 2. As moving 
u — > u n , one randomly chooses Ae = 1 or 2 and pro- 
pose to update the edge-ira n state as mod s(e + Ae). 
In other words, an 'empty' edge is proposed to randomly 
become a 'bond' or a 'dimer'; 'dimer' is to be 'empty' or 
'bond'; and 'bond' is to be 'empty' or 'dimer'. However, 
not all the proposals will generate a valid configuration 
that has at most two states in Fig. [fo] and the others in 
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Fig. [4] For instance, (1), in the M sector, when u = v 
is at state- 7 and the empty edge is proposed to become 
a dimer, the resulting vertex state at u will not be in 
Fig-IH (2), in the W sector, when u is at state-9' and the 
proposal is e = — > e = 1, this would yield state-11' at u 
which is not in Fig. 0] as required. A proposal would be 
rejected if it leads to an invalid configuration. On this 
basis, a version of the worm algorithm can be formulated 
as (Algorithm 2) 

1. If u = v, move it to a randomly chosen vertex. 

2. Same as in Algorithm 1. 

3. Same as in Algorithm 1. 

4. Randomly choose Ae = 1 or 2, and propose to up- 
date the edge-MM„ state as e UUn — > mod 3(e UUn + 
Ae). The proposal will be rejected if it yields 

• for u = v,V v or V Un {9', • • • , 13'} in Fig. El 

• for u ^ v and v ^ u n , V u ^ {6, 7} in Fig. @]or 
V Un £{9',-- - ,13'}; 

• for u ^= v and v — u n , V u or V Un ^ {6, 7}. 

In this case, step-5 will be skipped. Symbol V u 
represents the vertex state at u. 

5. Accept the update with probability 

V a = min [l,n A « (W^W^) / (W^W^})] , 

where A£ denotes the change of the loop number 
in the update. We remind that the constant £, w in 
Eq. ([11]) is set £ w = 1. 

Simulation consists of repetition of these steps, and the 
measurement is taken in the M sector. 

A practically important matter for implementing Al- 
gorithm 2 is that a non-local query is needed to calcu- 
late the loop-number difference A£. We shall follow the 
simultaneous breadth-first searching technique and the 
trick to avoid as much as possible queries, as described 
inRef. 

More importantly, one can apply the so-called color- 
ing method to avoid altogether the need for such global 
queries for n > 1. The key ingredient of the coloring 
method is the trivial identity n = 1 + (n — 1) for the 
statistical weight n of each loop. One can introduce an 
auxiliary variable c = 0, 1 and rewrite the identity as 



[18 cfi + (n - l)5 c>l ] 



(17) 



c=0,l 



The variable c is generally referred to as the coloring 
variable, and c = (1) is said 'active' ('inactive') . See 
Refs. [20] for details. In practise, the coloring variable 
is assigned to each vertex in the M sector as (Coloring 
assignment) 

1. Set all vertices off loops be active (c = 0). 



2. Independently for each loop, choose c = with 
probability 1/n and c = 1 with probability (1 — 
1/n), and assign it to all the vertices on the loop. 

On the basis of the Coloring assignment, the whole lattice 
G is divided into the active sublattice G a and the inac- 
tive sublattice Gi. In G a the vertices are active and the 
edges connect two active vertices; in Gi the vertices are 
inactive and the edges connect two inactive vertices. The 
edges connecting one active and one inactive vertex form 
the boundaries separating G a and Gi. We state that, 
conditioning on this decomposition, the vertex-state con- 
figuration on the induced sublattice G a and Gi is nothing 
but a generalized FPLD model with n' = 1 and (n — 1), 
respectively. 

One has now the right to update these generalized 
FPLD models via any valid Monte Carlo algorithm. We 
choose Algorithm 2 to update the model with n! — 1 on 
G a and the identity operation ('do nothing') on Gi. Due 
to the fact n' = 1, the loop- number change A£ does not 
matter anymore. Therefore, one can formulate another 
version of the worm algorithm as (Algorithm 3) 

1. Do the Coloring assignment if u = v. 

2. Do M times of the coarse Monte Carlo steps (from 
and back to the M sector) by performing Algorithm 
2 on the induced subgraph G a with n' = 1. 

The parameter M > 1 can be set such that step 1 and 2 
take comparable CPU time. 

For the actual implementation of Algorithm 2 and 3, 
positive statistical weights have to be assigned to vertex 
states in Fig. [6l Before discussing on this, we mention 
that there exist some freedom to choose which vertex 
state is allowed in the W sector. As long as ergodicity is 
satisfied, the consideration is to optimize the efficiency. 
In Fig. [6l we do not allow the state with two bonds and 
a dimer, because the only way to generate this state is 
to add a dimer to state- 7 and the only way to return to 
Fig. 0] is to delete the newly generated dimer. Thus, such 
a state will not help updating the configuration while 
increasing computational burden. In contrast, state-9' 
and -10' (-11' and -13') are important for moving around 
the dimers (bonds). We set 



W 9 , = Ww = l Wi2> = min(l, W 7 ) and 
Wxv = Wxv =W 7 =e 



-2K 



(18) 



State-12' is useful for switching between dimer and bond, 
but should not occur more frequently than state-6 or -7. 



IV. RESULTS 

The complete phase diagram of AT model on the tri- 
angular lattice is shown in Fig. [7l In following, we shall 
present numerical results and discuss the phase bound- 
ary in the antiferromagnetic two-spin coupling region 
(J<0). 
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FIG. 7: (Color online) Phase diagram of the AT model on 
the triangular lattice. Points-'P', T, 'O '-correspond to the 
4-state Potts, the Ising, and the 0(2) loop model, respec- 
tively. Points-'F', 'C, 'D'-denote the zero-temperature Ising 
antiferromagnet in variable ar, a or r (decoupled), a or r 
(correlated), respectively. Point 'E' is the zero-temperature 
4-state Potts antiferromagnet. 



A. Finite J 

We employ Algorithm 1 with Eq. (|14l) to simulate the 
AT model in the region of finite J < on triangular 
lattices with periodic boundary conditions, using system 
sizes in the range 6 < L < 192. 

For a given loop configuration, we generate the associ- 
ated spin configuration on the triangular lattice accord- 
ing to the low-temperature expansion rule. Note that, 
due to the periodic boundary condition, a loop config- 
uration may correspond to no spin configuration. This 
occurs when there exists an odd number of red or blue 
loops winding around the boundary. In this case, we take 
no measurement, and the simulation continues until the 
next try. Let X be the indicator function which is 1 if 
the loop configuration is measuring and corresponds to a 
spin configuration, and otherwise; let O be the operator 
computed in one of the 4 compatible spin configurations; 
therefore, what we are computing is [0AT]/[X], with [ | 
the statistical average over loop configurations. The non- 
valid loop configuration is not weighted, thus it does not 
influence any of the numerical data related to the spin 
variables. Further, since the special cases that the loop 
configuration docs not correspond to any spin configu- 
ration result from boundary effects, such cases do not 
dominate in large systems. 

Two types of magnetization are measured as 



and M a 



(19) 
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FIG. 8: (Color online) Quantity Q aT versus K at J = —1.0. 
Lines connecting the data points are for illustration purpose. 



ingly, the susceptibilities are defined as 

X* = V{Ml) and X *r = V(M, 



2 ) 

or I i 



(20) 



with ( ) for statistical average. Dimensionless ratios are 
found to be very powerful in locating the critical points 
of many systems under continuous phase transitions. On 
the basis of the fluctuation of the magnetization, we de- 
fine two distinct dimensionless ratios as 1251 



2\2 



and Q a 



(M, 



2 \2 
ari 



We also measure energy-like quantities as 

Ecr = -jy^ g t CTj 

(<,i> 

E = E a + E T + E aT , 



(21) 

(22) 

(23) 
(24) 



as well as the associated specific-heat-like quantities 
Ca = {{El) - {Ea) 2 )/V, Car, and C. 

The AT model for J = reduces to the standard Ising 
model in the Ising-spin variable err, and undergoes a 
Ising-like transition at K c . For K < K c , the configu- 
rations in the Ising variables a, r, and err are all in the 
disordered (paramagnetic) state; for K > K c , ar is in the 
ferromagnetic state while a and r are still in the para- 
magnetic state. We expect that this scenario continues 
into the region J < 0. 

We choose J = —0.2, —0.6, -1.0, and —2.0, and per- 
form some preliminary and coarse simulations to approx- 
imately locate the intersection of Q ar for various linear 
system sizes L. Then, fine and extensive simulations are 
carried out near the estimated critical point. Figure [8] 
displays Q aT versus K for different L at J = —1.0, indi- 
cating a critical point near K « 0.1265. 

The finite-size scaling behavior of Q aT (K, L) near the 
critical point K c is described by 



where the summation is over the whole lattice. Accord- 



Q(K,L) =Q(tL Vt ,bL Vi ) 



(25) 
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J 


-0.2 -0.6 -1.0 -2.0 


K e 


0.25303 (2) 0.18164(2) 0.12653 (2) 0.06306 (3) 




48 48 48 48 


xVdof 


1.05 0.86 1.14 1.21 



TABLE I: Details in the data fits according to Eq. [26] 



where t and i represent the leading and the subleading 
thermal scaling fields, with t oc (K — K c ) + ■ ■ ■ . The 
associated renormalization exponents are denoted as y t 
and yi. A Taylor expansion of Eq. ([23]) yields 26] 



Q(K,L) = Q c + ai AKL v * 
+ cAKL yt+y > +.. 



«2 



(AK) 2 L 2yt + bL Vi 



(26) 



with AK = K — K c . Parameters Oi, a 2 , b, and c are 
unknown constants. 

According to the least-squares criterion, we fit the Q aT 
data to Eq. (|2"6")l . Assuming the transition is Ising-likc, 
we expect that the leading two finite-size correction expo- 
nents are yi = 2-2y h = -7/4 and y 2 = yt — -2 for Q aT , 
where y^ = 15/8 is the magnetic renormalization expo- 
nent. With yi and y 2 fixed and L > L m j n = 48, we obtain 
K c = 0.12653(2), y t = 1.01(2), and Q c = 0.8587(1) for 
J = —1.0. The chi square per degree of freedom (x 2 /dof) 
is 1.14. The estimate of yt is consistent with the exact 
result yt = 1, and the universal ratio Q c — 0.8587 also 
agrees well with the earlier estimate Q c = 0.858 725 28(3) 
for the Ising model on the triangular lattice [27$ . 

The data of susceptibility x<tt is analyzed by 

X(K,L) = L- 2yh+d (a + ai AKL Vt + a 2 (AK) 2 L 2yt 
+ bL Vi + cAKL Vt+Vi + ...) , (27) 

and we determine the magnetic exponent as yh — 
1.876(2), in good agreement with the exact value yt = 
15/8. The specific-heat-like quantity C is also found to 
diverge approximately in the logarithmic scale as L in- 
creases. No phase transition is observed for Ising variable 
a or t. 

Similar results arc found for other values of J, and the 
estimated critical points are listed in Table |U 

On this basis, we conclude that the phase transition 
of the AT model in region (K > 0, J < 0) with finite J 
is in the Ising universality. Finally, we mention that the 
worm-type algorithm hereby does not suffer much from 
critical slowing-down. 



B. 



J — > — oo 



Table [T] suggests that the critical coupling K c becomes 
smaller as J becomes more negative, and that the ending 
point of the critical line for J — > — oo is very close to 
K = 0, since K C (J = -2) = 0.06306(3) is already near 0. 
To locate the ending point more accurately, we directly 
simulate the J — > — oo limit, which makes use of the 
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FIG. 9: (Color online) Dimensionless ratio Q aT versus K for 
the n=2 FPLD model. 



exact mapping to the n — 2 FPLD model and employs 
Algorithm 3. System sizes take 6 values in range 30 < 
L < 960. 

Note that the loops in the FPLD model serve as do- 
main walls for the Ising variable <jt in the AT model. 
According to the low-temperature expansion rule, on the 
triangular lattice we sample magnetization density M aT , 
susceptibility x<jt, dimensionless ratio Q aT , energy E aT , 
and specific heat C„t, whose definitions can be found in 
Eqs. (TT9")) - (|2"4"|) . Further, to explore the loop-length distri- 
bution, on the honeycomb lattice we measure the length 
of the longest loop as Si . 

The finite-size data of the dimensionless ratio Q aT are 
plotted in Fig. an eye- view fitting yields a critical point 
as K c = 0.00(2). For K > K c , the Q aT value rapidly ap- 
proaches to 1 as size L increases. This reflects that the 
Ising variable err exhibits a long-range ferromagnetic or- 
der on the triangular lattice; correspondingly, on the hon- 
eycomb lattice loops are small-i.e., in a disordered state. 
For K < K c , Q aT converges to a constant Q c which 
deviates from the trivial Gaussian value 1/3. This im- 
plies that, despite the absence of a long-range order, the 
spin-spin correlation function decays algebraically over 
the distance. 

In Fig. [HI one can observe that Q aT at K < rapidly 
converges to a K-dependent value, as expected in the low- 
temperature BKT phase. This reminds us the analogy 
between the FPLD and Nienhuis's 0(n) honeycomb loop 
model with n = 2. The phase diagram of the latter is 
shown in Fig. [TUl where x is the statistical weight for an 
occupied bond. For a given < n < 2, the O(n) loop 
model exhibits three distinct phases: a dilute and dis- 
ordered phase (small x), a densely-packed phase (large 
finite x), and a fully-packed phase (infinite x). Further- 
more, the model is exactly solvable on the curves [23[ 



1 

x± 



2± V2" 



(28) 



The system is equivalent to the tricritical q — n 2 Potts 
model along the critical line x + , belongs to the criti- 
cal q = n 2 Potts universality class in the densely-packed 



0.0 



n 



FIG. 10: (Color online) Phase diagram of the O(n) loop 
model [231] . Red lines denote the directions of the renormal- 
ization flows. 



phase, and is in another critical universality in the fully- 
packed phase. For n = 2, the two solvable lines x± 
merge at a single point; the renormalization field is 
marginally relevant (irrelevant) for x < x± (x > x±). In 
other words, the phase transition at x± is Berezinsky- 
Kosterlitz-Thouless(BKT)-like. At the special point 
x±(n — 2), the amplitude of the renormalization field is 
zero, and thus logarithmic corrections, present at most of 
BKT-like critical points, disappear. This explains the ab- 
sence of logarithmic corrections in the critical Baxter- Wu 
model, which can be exactly mapped onto the 0(2) loop 
model at x± . For the critical 0(2) loop model, it has been 
identified that Si oc L yH = L 3 / 2 and X *r oc L 2ym ~ 2 = L, 
where yu = 3/2 is the hull exponent and yto = 3/2 is 
the leading thermal renormalization exponent in the lan- 
guage of the Potts model [28j . 

Since the state space of the FPLD model is a subspace 
of the 0(2) loop model, it is reasonable to conjecture 
that the two models are in the same universality class. 
Namely, we expect that the FPLD model undergoes a 
BKT-like transition at K c , where the logarithmic correc- 
tions are absent; for K < K c the system is in the same 
universality class as at K c but with logarithmic correc- 
tions; for K — > — oo it is in another universality class. 
Making use of the known exponent yn = 3/2 for Si and 
2y t a — 2 = 1 for Xar, we plot L~ 3 / 2 S\ and L~ 1 \ rjT versus 
K in Figs. QT] and [12] respectively. They both display 
a nice intersection at K = 0.000. From Fig. [T^Jone can 
observe that the exponent yn varies along the BKT crit- 
ical line, which reconciles the difference of yn between 
the present model and the two-dimensional XY models. 



To further explore the potential logarithmic correc- 
tions, we assume K c = and plot L~ 3 l 2 Si and L~ 1 x<jt 
at K — versus L^ 1 . As shown in Fig.[T31 the rapid con- 
vergence implies the absence of logarithmic corrections; 
corrections with term L _1 are also very weak if they ex- 
ist. 

According to the least-squares criterion, the Si and 
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FIG. 11: (Color online) L~ 3/2 S! versus K for the n=2 FPLD 
model. 
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FIG. 12: (Color online) IT x Xtt versus K for the n=2 FPLD 
model. 



Xctt data are fitted by 

Y(K, L) = c + ci AK + ■ ■ ■ + L Xy (a + ai AK In L 
+a 2 (AK) 2 In 2 L + b r L yi + b 2 L y * + ...) . (29) 

Here a, are coefficients of the finite-size scaling variable 
AK In L with AK = K — K c , hi are amplitudes of finitc- 
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FIG. 13: (Color online) Quantities L~ 3 ^ 2 Si — ao, sl and 
L~ 1 Xot — ao,x at A" = versus L _1 for the n=2 FPLD model. 
Constants ao, 31 and ao iX are obtained from the fits. 
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size corrections, and the terms with a accounts for ana- 
lytical background. There are also cross-terms involving 
products of terms arising from these three sources. The 
exponent Xy is a general label for quantity Y. Equa- 
tion ([29]) has assumed the absence of logarithmic cor- 
rections. It occurs that both the Si and the x^t data 
with L > L m ; n = 60 can be well described by Eq. ([25)1 
with hxed correction exponents yi = — 1 and j/2 = —2. 
The results for Si are Xy = 1.498(3), K c = 0.001(2), 
and x 2 /dof=1.22; for X ar are X Y = 1.001(2), K c = 
-0.001(1), and x 2 /dof=0.87. These agree well with the 
known exponents 3/2 for S\ and 1 for Xar , as well as with 
the expectation K c — 0. If the exponents Xy are further 
fixed at the known values, we obtain K c = 0.0002(3), 
X 2 /dof=0.94 from Si and -0.0003(4), x 2 /dof=1.09 from 
Xar- On this basis, we estimate the critical point as 
A c = —0.0001(6), which covers the uncertainties of K c 
from Si and Xa-r- 

We mention that, when an external field of strength 
h/T is applied to the triangular Ising antiferromagnet, 
the critical state of the system is not immediately de- 
stroyed. Instead, the system has a BKT-like transition 
at h c = 0.266(10) (29j. However, our Monte Carlo results 
suggest that a critical point A c ^ is rather unlikely for 
the n = 2 FPLD model. 

In the limit A — > oo, the Ising variable err is in the 
ferromagnetic state. However, in terms of the a or the 
r variable, it can be easily derived that the system is 
also an Ising model with coupling 2 J. Namely, along the 
tanh A' = 1 line, the AT model has an Ising-like tran- 
sition at tanh2J = y2 — 1. Further, the corner point 
D := (tanhA = l,tanhJ = —1) corresponds to the 
triangular antiferromagnet at zero temperature, which 
is critical. Together with the earlier discussions in Sec. 
I, this means that, in Fig. the limiting points-D, C, 
O, F, E-are all critical. From our simulations in range 
-0.2 < K < 0.1 along the tanh J = -1 line (EC+CD), 
we observe that, in the whole range, there exist alge- 
braically decaying two-point correlation function for the 
a or the r variable. On this basis, we conjecture that the 
whole tanh J = — 1 line (EC+CD) is critical for the a or 
the t variable. Simulation along the tanh A = — 1 line us- 
ing the present worm algorithms suffers significantly from 
critical slowing-down. Nevertheless, we suspect that the 
whole tanh A = — 1 line is critical for the o~t variable. 



V. DYNAMIC CRITICAL BEHAVIOR 

In this section, we briefly report the efficiency of Algo- 
rithm 2 for the n = 2 FPLD model, using the standard 
procedure described in Ref. [30j. 

For each observable (say O), we calculate its autocor- 
relation function 

p o (t) = (O(t)O(0))-(Or, 

where ( ) denotes expectation with respect to the sta- 
tionary distribution. We then obtain the corresponding 




FIG. 14: (Color online) Ln(T B /L 2 ) versus InL at K = 0. 



integrated autocorrelation time as 



(30) 



The dynamic critical exponent z- lnt Q is defined by 

T in t,o - r nt -° • (31) 

where £ is the spatial correlation length. On a finite lat- 
tice at criticality, £ is cut off by system size L. Therefore, 
one has 



r mt,e> 



bL z 



(32) 



with a and b unknown parameters. 

We simulate at the critical point A c = 0. Note that, 
during the worm simulations we measure the observables 
only when the chain visits the Eulerian subspace, roughly 
every Tg ~ T J d-2X c However, it is natural to de- 

fine Zj n t,0 as in Ref. [l9| to measure time in units of 
sweeps of the lattice, i.e. L d hits. Since one sweep takes 
of order L 2X " visits to the Eulerian subspace, we have 
t ~ L z+2Xc . As shown in FigfLlI the exponent 2A e is 
estimated to be 0.50(1). 

Among the measured quantities including the longest- 
loop length Si, the loop number t, and the energy- like 
quantity E aT etc, E aT is found to have the largest value 
of r int . Figure displays PE aT {t/T int ^ EaT ) as a func- 
tion of t/Ti n t,E aT , where an approximately exponential 
decay is observed. The Ti nti E aT data are analyzed, and 
we obtain z- ln t,E aT = 0.28(1), which is shown in Fig. [T^l 
Similar fits are done for other quantities, and we have 
z int Sl = 0.26(1) and z int j = 0.27(1). In these fits, x 2 /dof 
ranges from 0.74 to 1.31. Therefore, our numerical results 
suggest that the present worm algorithm is even more ef- 
ficient than the one in Ref. [19|. 

Simulations are also carried out for A = —0.05, and 
the dynamic critical behavior cannot be distinguished 
from that for A = 0. 
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FIG. 15: (Color online) pE < , T {t/T int ,E aT ) versus t/T int .E aT at 
K = 0. 
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FIG. 16: (Color online) Ln(ri n t,o) versus InL for different 
observables at K — 0. 



VI. DISCUSSION 

In summary, we have formulated two versions of the 
worm-type algorithms for the AT model on the triangular 
lattice. The algorithms are based on the low-temperature 
expansion graph of the AT model, and use the language of 
vertex states. Such a formulation not only provides us a 
different angle to understand the worm method, but also 
offers an easy way to overcome the absorbing difficulty. 
The efficiency of our algorithm is studied and can also 
be reflected by the fact that we can simulate up to size 
L = 960. Apparently, Algorithm 1 can be applied to the 
ferromagnetic region J > of the triangular AT model 
and to the AT model on other planar lattices. Further, we 
mention that the worm-type algorithms can be developed 
on the basis of the high-temperature expansion graph of 
the AT model. This yields a graphical model also by 
Eq. ([3]), but defined on the original lattice for the AT 
model. The statistical weights of the occupied bonds are 

X r = X b = (e 2X sinh2J)/(e 2 * r cosh2J + 1) 
X r+b = (e 2i ^cosh2J- l)/(e 2K cosh2J + 1) . (33) 

It is reasonable to expect good efficiency for the AT 
model on non-planar lattices-e.g., in higher spatial 



J 


-0.2 -0.4 -0.6 -1.0 -2.0 -oo 




0.442(2) 0.408(2) 0.386(2) 0.370(2) 0.366(2) 0.3655(3) 


Lmin 


36 36 36 36 36 36 


xVdof 


1.21 0.91 1.07 1.23 0.85 1.15 



TABLE II: Details in the data fits according to Eq.[26]on the 
kagome lattice. 
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FIG. 17: (Color online) Phase diagram of the AT model on 
the kagome lattice. 



dimensions-with non-negative weights X r = X\, and 

X r +b- 

The high efficiency of the worm algorithms allows us 
to explore the triangular-lattice AT model in the antifcr- 
romagnetic region, and accordingly we conjecture a com- 
plete phase diagram in the (J, K) plane. Of the particu- 
lar interest is the J — > — oo limit, where the AT model is 
mapped onto the FPLD model with n = 2. As suggested 
by the Monte Carlo simulation, the AT model undergoes 
a BKT-like transition along the tanh J = — 1 line, in the 
same universality class as the classical XY model. We 
also mention that it remains to be explored whether or 
not, for other values of n, the FPLD and Nienhuis's O(n) 
model are in the same universality class. 

Finally, we perform simulations for the AT model on 
the kagome lattice in the region (J < Q,K > 0), and 
determine a line of Ising-like critical points. The results 
are shown in Table [TTJ Unlike on the triangular lattice, 
the critical line ends at K c = 0.3655 > 0, still in the Ising 
universality. Taking into account that the frustration on 
the kagome lattice is only partial, this is not surprising. 
Accordingly, the phase diagram is shown in Fig. 1171 
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